function TL = simu_timelines(n0, n, T, c)
% Simulate timelines
%
%   TL = simu_timelines(n0, n, T, c);
%       simulate timelines (a positive value that evolve over time).
%
%       Inputs:
%       - n0:       the number of initially existent timelines
%       - n:        the number of timelines to generate on the fly
%       - T:        the length of observed period
%       - c:        the configuration control the span of each timeline.
%                   It is a triple in the following format:
%
%                       (tt, bt, vt)
%                   
%                   tt:  the transit time
%                   bt:  the base time (fixed part)
%                   vt:  the average length of the variable part
%
%       The output TL is an (n+n0) x (T+1) matrix representing the 
%       timelines. In particular, each row of TL is a timeline, each 
%       column of TL gives the values of different timelines at a 
%       particular time.
%

%   Created by Dahua Lin, on Nov 24, 2010
%

%% verify input argument

if ~(isnumeric(n0) && isscalar(n0) && n0 == fix(n0) && n0 >= 0)
    error('simu_timelines:invalidarg', ...
        'n0 should be a non-negative integer scalar.');
end

if ~(isnumeric(n) && isscalar(n) && n == fix(n) && n >= 0)
    error('simu_timelines:invalidarg', ...
        'n should be a non-negative integer scalar.');
end

if ~(isnumeric(T) && isscalar(T) && T == fix(T) && T >= 1)
    error('simu_timelines:invalidarg', ...
        'T should be a positive integer scalar.');
end

if ~(isnumeric(c) && numel(c) == 3)
    error('simu_timelines:invalidarg', 'c should be a numeric triple.');
end
   

%% main

TL = zeros(n0+n, T+1);

for i = 1 : n0
    TL(i, :) = make_tl(0, T, c);
end

for i = 1 : n
    TL(n0+i, :) = make_tl(rand() * T, T, c);
end


%% core sub-routines

function tl = make_tl(s, T, c)

tt = c(1); bt = c(2); vt = c(3);

tl = zeros(1, T+1);

% fade-in
if s > 0 && s < T
    t = ceil(s) : min(floor(s+tt), T);
    if ~isempty(t)
        tl(t+1) = (t - s) / tt;
    end
    s = s + tt;
end

% main duration

if s < T
    dt = bt - log(rand()) * vt;
    t = ceil(s) : min(floor(s+dt), T);
    if ~isempty(t)
        tl(t+1) = 1;
    end
    s = s + dt;
end

% fade-out
if s < T
    t = ceil(s) : min(floor(s+tt), T);
    if ~isempty(t)
        tl(t+1) = (s+tt-t) / tt;
    end
end

